//
// File: freqDomainHRV.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 23-Apr-2025 16:26:18
//

// Include Files
#include "freqDomainHRV.h"
#include "blockedSummation.h"
#include "colon.h"
#include "combineVectorElements.h"
#include "detrend.h"
#include "freqDomainHRV_data.h"
#include "freqDomainHRV_initialize.h"
#include "freqDomainHRV_rtwutil.h"
#include "freqDomainHRV_types.h"
#include "hamming.h"
#include "interp1.h"
#include "linspace.h"
#include "lomb2.h"
#include "minOrMax.h"
#include "pburg.h"
#include "pwelch.h"
#include "rt_nonfinite.h"
#include "trapz.h"
#include "coder_array.h"
#include <cmath>
#include <string.h>

// Function Declarations
static void b_times(coder::array<double, 1U> &in1,
                    const coder::array<double, 1U> &in2);

static void calcAreas(const coder::array<double, 1U> &F,
                      const coder::array<double, 2U> &PSD, const double VLF[2],
                      const double LF[2], const double HF[2],
                      struct2_T *output);

static void calcAreas(const coder::array<double, 2U> &F,
                      coder::array<double, 2U> &PSD, const double VLF[2],
                      const double LF[2], const double HF[2],
                      struct2_T *output);

// Function Definitions
//
// Arguments    : coder::array<double, 1U> &in1
//                const coder::array<double, 1U> &in2
// Return Type  : void
//
static void b_times(coder::array<double, 1U> &in1,
                    const coder::array<double, 1U> &in2)
{
  coder::array<double, 1U> b_in1;
  int i;
  int loop_ub;
  int stride_0_0;
  int stride_1_0;
  if (in2.size(0) == 1) {
    i = in1.size(0);
  } else {
    i = in2.size(0);
  }
  b_in1.set_size(i);
  stride_0_0 = (in1.size(0) != 1);
  stride_1_0 = (in2.size(0) != 1);
  if (in2.size(0) == 1) {
    loop_ub = in1.size(0);
  } else {
    loop_ub = in2.size(0);
  }
  for (i = 0; i < loop_ub; i++) {
    b_in1[i] = in1[i * stride_0_0] * in2[i * stride_1_0];
  }
  in1.set_size(b_in1.size(0));
  loop_ub = b_in1.size(0);
  for (i = 0; i < loop_ub; i++) {
    in1[i] = b_in1[i];
  }
}

//
// calcAreas - Calulates areas/energy under the PSD curve within the freq
// bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
// percentage, and normalized units. Also returns LF/HF ratio.
//
// Inputs:
//    PSD: PSD vector
//    F: Freq vector
//    VLF, LF, HF: array containing VLF, LF, and HF freq limits
//    flagNormalize: option to normalize PSD to max(PSD)
// Output:
//
// Usage:
//
//
//    Modified from Gary Clifford's ECG Toolbox: calc_lfhf.m
//
// Arguments    : const coder::array<double, 1U> &F
//                const coder::array<double, 2U> &PSD
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                struct2_T *output
// Return Type  : void
//
static void calcAreas(const coder::array<double, 1U> &F,
                      const coder::array<double, 2U> &PSD, const double VLF[2],
                      const double LF[2], const double HF[2], struct2_T *output)
{
  coder::array<double, 1U> b_F;
  coder::array<double, 1U> b_PSD;
  coder::array<double, 1U> tmpF;
  coder::array<int, 1U> r10;
  coder::array<int, 1U> r11;
  coder::array<int, 1U> r6;
  coder::array<int, 1U> r7;
  coder::array<int, 1U> r8;
  coder::array<int, 1U> r9;
  coder::array<boolean_T, 1U> r;
  coder::array<boolean_T, 1U> r1;
  coder::array<boolean_T, 1U> r2;
  coder::array<boolean_T, 1U> r3;
  coder::array<boolean_T, 1U> r4;
  coder::array<boolean_T, 1U> r5;
  double aLF;
  double aVLF;
  double peakLF;
  double peakVLF;
  int end;
  int iindx;
  int loop_ub;
  // normalize PSD if needed
  //  find the indexes corresponding to the VLF, LF, and HF bands
  r.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r[end] = (F[end] >= VLF[0]);
  }
  r1.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r1[end] = (F[end] <= VLF[1]);
  }
  r2.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r2[end] = (F[end] >= LF[0]);
  }
  r3.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r3[end] = (F[end] <= LF[1]);
  }
  r4.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r4[end] = (F[end] >= HF[0]);
  }
  r5.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r5[end] = (F[end] <= HF[1]);
  }
  // Find peaks
  // VLF Peak
  end = r.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  tmpF.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      tmpF[loop_ub] = F[i];
      loop_ub++;
    }
  }
  end = r.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  r6.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      r6[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  b_PSD.set_size(r6.size(0));
  loop_ub = r6.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r6[end] - 1];
  }
  coder::internal::maximum(b_PSD, &aVLF, &iindx);
  peakVLF = tmpF[iindx - 1];
  // LF Peak
  end = r2.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  tmpF.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      tmpF[loop_ub] = F[i];
      loop_ub++;
    }
  }
  end = r2.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  r7.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      r7[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  b_PSD.set_size(r7.size(0));
  loop_ub = r7.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r7[end] - 1];
  }
  coder::internal::maximum(b_PSD, &aVLF, &iindx);
  peakLF = tmpF[iindx - 1];
  // HF Peak
  end = r4.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  tmpF.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      tmpF[loop_ub] = F[i];
      loop_ub++;
    }
  }
  end = r4.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  r8.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      r8[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  b_PSD.set_size(r8.size(0));
  loop_ub = r8.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r8[end] - 1];
  }
  coder::internal::maximum(b_PSD, &aVLF, &iindx);
  //  calculate raw areas (power under curve), within the freq bands (ms^2)
  end = r.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  r9.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      r9[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  b_F.set_size(r9.size(0));
  loop_ub = r9.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_F[end] = F[r9[end] - 1];
  }
  b_PSD.set_size(r9.size(0));
  loop_ub = r9.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r9[end] - 1];
  }
  aVLF = coder::trapz(b_F, b_PSD);
  end = r2.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  r10.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      r10[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  b_F.set_size(r10.size(0));
  loop_ub = r10.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_F[end] = F[r10[end] - 1];
  }
  b_PSD.set_size(r10.size(0));
  loop_ub = r10.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r10[end] - 1];
  }
  aLF = coder::trapz(b_F, b_PSD);
  end = r4.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  r11.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      r11[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  b_F.set_size(r11.size(0));
  loop_ub = r11.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_F[end] = F[r11[end] - 1];
  }
  b_PSD.set_size(r11.size(0));
  loop_ub = r11.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r11[end] - 1];
  }
  double aHF;
  double aTotal;
  aHF = coder::trapz(b_F, b_PSD);
  // aVLF=aVLF/1e12;
  // aLF=aLF/1e12;
  // aHF=aHF/1e12;
  aTotal = (aVLF + aLF) + aHF;
  // calculate areas relative to the total area (%)
  // calculate normalized areas (relative to HF+LF, n.u.)
  // calculate LF/HF ratio
  // create output structure
  output->aVLF = std::round(aVLF * 100.0) / 100.0;
  //  round
  output->aLF = std::round(aLF * 100.0) / 100.0;
  output->aHF = std::round(aHF * 100.0) / 100.0;
  output->aTotal = std::round(aTotal * 100.0) / 100.0;
  output->pVLF = std::round(aVLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pLF = std::round(aLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pHF = std::round(aHF / aTotal * 100.0 * 10.0) / 10.0;
  aVLF = aLF + aHF;
  output->nLF = std::round(aLF / aVLF * 1000.0) / 1000.0;
  output->nHF = std::round(aHF / aVLF * 1000.0) / 1000.0;
  output->LFHF = std::round(aLF / aHF * 1000.0) / 1000.0;
  output->peakVLF = std::round(peakVLF * 100.0) / 100.0;
  output->peakLF = std::round(peakLF * 100.0) / 100.0;
  output->peakHF = std::round(tmpF[iindx - 1] * 100.0) / 100.0;
}

//
// calcAreas - Calulates areas/energy under the PSD curve within the freq
// bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
// percentage, and normalized units. Also returns LF/HF ratio.
//
// Inputs:
//    PSD: PSD vector
//    F: Freq vector
//    VLF, LF, HF: array containing VLF, LF, and HF freq limits
//    flagNormalize: option to normalize PSD to max(PSD)
// Output:
//
// Usage:
//
//
//    Modified from Gary Clifford's ECG Toolbox: calc_lfhf.m
//
// Arguments    : const coder::array<double, 2U> &F
//                coder::array<double, 2U> &PSD
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                struct2_T *output
// Return Type  : void
//
static void calcAreas(const coder::array<double, 2U> &F,
                      coder::array<double, 2U> &PSD, const double VLF[2],
                      const double LF[2], const double HF[2], struct2_T *output)
{
  coder::array<double, 2U> b_PSD;
  coder::array<double, 2U> tmpF;
  coder::array<double, 2U> tmppsd;
  coder::array<int, 2U> r6;
  coder::array<int, 2U> r7;
  coder::array<int, 2U> r8;
  coder::array<boolean_T, 2U> r;
  coder::array<boolean_T, 2U> r1;
  coder::array<boolean_T, 2U> r2;
  coder::array<boolean_T, 2U> r3;
  coder::array<boolean_T, 2U> r4;
  coder::array<boolean_T, 2U> r5;
  double aLF;
  double aTotal;
  double aVLF;
  double peakLF;
  double peakVLF;
  int i;
  int idx;
  int iindx;
  int k;
  int last;
  // normalize PSD if needed
  last = PSD.size(1);
  if (PSD.size(1) <= 2) {
    if (PSD.size(1) == 1) {
      aVLF = PSD[0];
    } else if ((PSD[0] < PSD[PSD.size(1) - 1]) ||
               (std::isnan(PSD[0]) && (!std::isnan(PSD[PSD.size(1) - 1])))) {
      aVLF = PSD[PSD.size(1) - 1];
    } else {
      aVLF = PSD[0];
    }
  } else {
    if (!std::isnan(PSD[0])) {
      idx = 1;
    } else {
      boolean_T exitg1;
      idx = 0;
      k = 2;
      exitg1 = false;
      while ((!exitg1) && (k <= last)) {
        if (!std::isnan(PSD[k - 1])) {
          idx = k;
          exitg1 = true;
        } else {
          k++;
        }
      }
    }
    if (idx == 0) {
      aVLF = PSD[0];
    } else {
      aVLF = PSD[idx - 1];
      i = idx + 1;
      for (k = i; k <= last; k++) {
        aTotal = PSD[k - 1];
        if (aVLF < aTotal) {
          aVLF = aTotal;
        }
      }
    }
  }
  PSD.set_size(1, PSD.size(1));
  idx = PSD.size(1) - 1;
  for (i = 0; i <= idx; i++) {
    PSD[i] = PSD[i] / aVLF;
  }
  //  find the indexes corresponding to the VLF, LF, and HF bands
  r.set_size(1, F.size(1));
  idx = F.size(1);
  for (i = 0; i < idx; i++) {
    r[i] = (F[i] >= VLF[0]);
  }
  r1.set_size(1, F.size(1));
  idx = F.size(1);
  for (i = 0; i < idx; i++) {
    r1[i] = (F[i] <= VLF[1]);
  }
  r2.set_size(1, F.size(1));
  idx = F.size(1);
  for (i = 0; i < idx; i++) {
    r2[i] = (F[i] >= LF[0]);
  }
  r3.set_size(1, F.size(1));
  idx = F.size(1);
  for (i = 0; i < idx; i++) {
    r3[i] = (F[i] <= LF[1]);
  }
  r4.set_size(1, F.size(1));
  idx = F.size(1);
  for (i = 0; i < idx; i++) {
    r4[i] = (F[i] >= HF[0]);
  }
  r5.set_size(1, F.size(1));
  idx = F.size(1);
  for (i = 0; i < idx; i++) {
    r5[i] = (F[i] <= HF[1]);
  }
  // Find peaks
  // VLF Peak
  last = r.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r[k] && r1[k]) {
      idx++;
    }
  }
  tmpF.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r[k] && r1[k]) {
      tmpF[idx] = F[k];
      idx++;
    }
  }
  last = r.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r[k] && r1[k]) {
      idx++;
    }
  }
  tmppsd.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r[k] && r1[k]) {
      tmppsd[idx] = PSD[k];
      idx++;
    }
  }
  coder::internal::maximum(tmppsd, &aVLF, &iindx);
  peakVLF = tmpF[iindx - 1];
  // LF Peak
  last = r2.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r2[k] && r3[k]) {
      idx++;
    }
  }
  tmpF.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r2[k] && r3[k]) {
      tmpF[idx] = F[k];
      idx++;
    }
  }
  last = r2.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r2[k] && r3[k]) {
      idx++;
    }
  }
  tmppsd.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r2[k] && r3[k]) {
      tmppsd[idx] = PSD[k];
      idx++;
    }
  }
  coder::internal::maximum(tmppsd, &aVLF, &iindx);
  peakLF = tmpF[iindx - 1];
  // HF Peak
  last = r4.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r4[k] && r5[k]) {
      idx++;
    }
  }
  tmpF.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r4[k] && r5[k]) {
      tmpF[idx] = F[k];
      idx++;
    }
  }
  last = r4.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r4[k] && r5[k]) {
      idx++;
    }
  }
  tmppsd.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r4[k] && r5[k]) {
      tmppsd[idx] = PSD[k];
      idx++;
    }
  }
  coder::internal::maximum(tmppsd, &aVLF, &iindx);
  //  calculate raw areas (power under curve), within the freq bands (ms^2)
  last = r.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r[k] && r1[k]) {
      idx++;
    }
  }
  r6.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r[k] && r1[k]) {
      r6[idx] = k + 1;
      idx++;
    }
  }
  tmppsd.set_size(1, r6.size(1));
  idx = r6.size(1);
  for (i = 0; i < idx; i++) {
    tmppsd[i] = F[r6[i] - 1];
  }
  b_PSD.set_size(1, r6.size(1));
  idx = r6.size(1);
  for (i = 0; i < idx; i++) {
    b_PSD[i] = PSD[r6[i] - 1];
  }
  aVLF = coder::trapz(tmppsd, b_PSD);
  last = r2.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r2[k] && r3[k]) {
      idx++;
    }
  }
  r7.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r2[k] && r3[k]) {
      r7[idx] = k + 1;
      idx++;
    }
  }
  tmppsd.set_size(1, r7.size(1));
  idx = r7.size(1);
  for (i = 0; i < idx; i++) {
    tmppsd[i] = F[r7[i] - 1];
  }
  b_PSD.set_size(1, r7.size(1));
  idx = r7.size(1);
  for (i = 0; i < idx; i++) {
    b_PSD[i] = PSD[r7[i] - 1];
  }
  aLF = coder::trapz(tmppsd, b_PSD);
  last = r4.size(1) - 1;
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r4[k] && r5[k]) {
      idx++;
    }
  }
  r8.set_size(1, idx);
  idx = 0;
  for (k = 0; k <= last; k++) {
    if (r4[k] && r5[k]) {
      r8[idx] = k + 1;
      idx++;
    }
  }
  tmppsd.set_size(1, r8.size(1));
  idx = r8.size(1);
  for (i = 0; i < idx; i++) {
    tmppsd[i] = F[r8[i] - 1];
  }
  b_PSD.set_size(1, r8.size(1));
  idx = r8.size(1);
  for (i = 0; i < idx; i++) {
    b_PSD[i] = PSD[r8[i] - 1];
  }
  double aHF;
  aHF = coder::trapz(tmppsd, b_PSD);
  // aVLF=aVLF/1e12;
  // aLF=aLF/1e12;
  // aHF=aHF/1e12;
  aTotal = (aVLF + aLF) + aHF;
  // calculate areas relative to the total area (%)
  // calculate normalized areas (relative to HF+LF, n.u.)
  // calculate LF/HF ratio
  // create output structure
  output->aVLF = std::round(aVLF * 1000.0) / 1000.0;
  output->aLF = std::round(aLF * 1000.0) / 1000.0;
  output->aHF = std::round(aHF * 1000.0) / 1000.0;
  output->aTotal = std::round(aTotal * 1000.0) / 1000.0;
  output->pVLF = std::round(aVLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pLF = std::round(aLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pHF = std::round(aHF / aTotal * 100.0 * 10.0) / 10.0;
  aTotal = aLF + aHF;
  output->nLF = std::round(aLF / aTotal * 1000.0) / 1000.0;
  output->nHF = std::round(aHF / aTotal * 1000.0) / 1000.0;
  output->LFHF = std::round(aLF / aHF * 1000.0) / 1000.0;
  output->peakVLF = std::round(peakVLF * 100.0) / 100.0;
  output->peakLF = std::round(peakLF * 100.0) / 100.0;
  output->peakHF = std::round(tmpF[iindx - 1] * 100.0) / 100.0;
}

//
// calcAreas - Calulates areas/energy under the PSD curve within the freq
// bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
// percentage, and normalized units. Also returns LF/HF ratio.
//
// Inputs:
//    PSD: PSD vector
//    F: Freq vector
//    VLF, LF, HF: array containing VLF, LF, and HF freq limits
//    flagNormalize: option to normalize PSD to max(PSD)
// Output:
//
// Usage:
//
//
//    Modified from Gary Clifford's ECG Toolbox: calc_lfhf.m
//
// Arguments    : const coder::array<double, 1U> &F
//                const coder::array<double, 1U> &PSD
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                struct2_T *output
// Return Type  : void
//
void calcAreas(const coder::array<double, 1U> &F,
               const coder::array<double, 1U> &PSD, const double VLF[2],
               const double LF[2], const double HF[2], struct2_T *output)
{
  coder::array<double, 1U> b_PSD;
  coder::array<double, 1U> tmpF;
  coder::array<double, 1U> tmppsd;
  coder::array<int, 1U> r6;
  coder::array<int, 1U> r7;
  coder::array<int, 1U> r8;
  coder::array<boolean_T, 1U> r;
  coder::array<boolean_T, 1U> r1;
  coder::array<boolean_T, 1U> r2;
  coder::array<boolean_T, 1U> r3;
  coder::array<boolean_T, 1U> r4;
  coder::array<boolean_T, 1U> r5;
  double aLF;
  double aVLF;
  double peakLF;
  double peakVLF;
  int end;
  int iindx;
  int loop_ub;
  // normalize PSD if needed
  //  find the indexes corresponding to the VLF, LF, and HF bands
  r.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r[end] = (F[end] >= VLF[0]);
  }
  r1.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r1[end] = (F[end] <= VLF[1]);
  }
  r2.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r2[end] = (F[end] >= LF[0]);
  }
  r3.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r3[end] = (F[end] <= LF[1]);
  }
  r4.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r4[end] = (F[end] >= HF[0]);
  }
  r5.set_size(F.size(0));
  loop_ub = F.size(0);
  for (end = 0; end < loop_ub; end++) {
    r5[end] = (F[end] <= HF[1]);
  }
  // Find peaks
  // VLF Peak
  end = r.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  tmpF.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      tmpF[loop_ub] = F[i];
      loop_ub++;
    }
  }
  end = r.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  tmppsd.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      tmppsd[loop_ub] = PSD[i];
      loop_ub++;
    }
  }
  coder::internal::maximum(tmppsd, &aVLF, &iindx);
  peakVLF = tmpF[iindx - 1];
  // LF Peak
  end = r2.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  tmpF.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      tmpF[loop_ub] = F[i];
      loop_ub++;
    }
  }
  end = r2.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  tmppsd.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      tmppsd[loop_ub] = PSD[i];
      loop_ub++;
    }
  }
  coder::internal::maximum(tmppsd, &aVLF, &iindx);
  peakLF = tmpF[iindx - 1];
  // HF Peak
  end = r4.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  tmpF.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      tmpF[loop_ub] = F[i];
      loop_ub++;
    }
  }
  end = r4.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  tmppsd.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      tmppsd[loop_ub] = PSD[i];
      loop_ub++;
    }
  }
  coder::internal::maximum(tmppsd, &aVLF, &iindx);
  //  calculate raw areas (power under curve), within the freq bands (ms^2)
  end = r.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  r6.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      r6[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  tmppsd.set_size(r6.size(0));
  loop_ub = r6.size(0);
  for (end = 0; end < loop_ub; end++) {
    tmppsd[end] = F[r6[end] - 1];
  }
  b_PSD.set_size(r6.size(0));
  loop_ub = r6.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r6[end] - 1];
  }
  aVLF = coder::trapz(tmppsd, b_PSD);
  end = r2.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  r7.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      r7[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  tmppsd.set_size(r7.size(0));
  loop_ub = r7.size(0);
  for (end = 0; end < loop_ub; end++) {
    tmppsd[end] = F[r7[end] - 1];
  }
  b_PSD.set_size(r7.size(0));
  loop_ub = r7.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r7[end] - 1];
  }
  aLF = coder::trapz(tmppsd, b_PSD);
  end = r4.size(0) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  r8.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      r8[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  tmppsd.set_size(r8.size(0));
  loop_ub = r8.size(0);
  for (end = 0; end < loop_ub; end++) {
    tmppsd[end] = F[r8[end] - 1];
  }
  b_PSD.set_size(r8.size(0));
  loop_ub = r8.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r8[end] - 1];
  }
  double aHF;
  double aTotal;
  aHF = coder::trapz(tmppsd, b_PSD);
  // aVLF=aVLF/1e12;
  // aLF=aLF/1e12;
  // aHF=aHF/1e12;
  aTotal = (aVLF + aLF) + aHF;
  // calculate areas relative to the total area (%)
  // calculate normalized areas (relative to HF+LF, n.u.)
  // calculate LF/HF ratio
  // create output structure
  output->aVLF = std::round(aVLF * 100.0) / 100.0;
  //  round
  output->aLF = std::round(aLF * 100.0) / 100.0;
  output->aHF = std::round(aHF * 100.0) / 100.0;
  output->aTotal = std::round(aTotal * 100.0) / 100.0;
  output->pVLF = std::round(aVLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pLF = std::round(aLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pHF = std::round(aHF / aTotal * 100.0 * 10.0) / 10.0;
  aVLF = aLF + aHF;
  output->nLF = std::round(aLF / aVLF * 1000.0) / 1000.0;
  output->nHF = std::round(aHF / aVLF * 1000.0) / 1000.0;
  output->LFHF = std::round(aLF / aHF * 1000.0) / 1000.0;
  output->peakVLF = std::round(peakVLF * 100.0) / 100.0;
  output->peakLF = std::round(peakLF * 100.0) / 100.0;
  output->peakHF = std::round(tmpF[iindx - 1] * 100.0) / 100.0;
}

//
// freqDomainHRV - calculates freq domain HRV using FFT, AR, and Lomb-Scargle
// methods
//
// Inputs:    ibi = 2Dim array of time (s) and inter-beat interval (s)
//            AR_order = order of AR model
//            window = # of samples in window
//            noverlap = # of samples to overlap
//            fs = cubic spline interpolation rate / resample rate (Hz)
//            nfft = # of points in the frequency axis
//            methods = methods of calculating freqDomain. Default is all 3
//                methhods
// Outputs:   output is a structure containg all HRV. One field for each PSD
// method
//            Output units include:
//                peakHF,LF,VLF (Hz)
//                aHF,aLF,aVLF (ms^2)
//                pHF,pLF,pVLF (%)
//                nHF,nLF,nVLF (%)
//                PSD (ms^2/Hz)
//                F (Hz)
// Usage:
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//  Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
//
//  This file is part of HRVAS
//
//  HRVAS is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  HRVAS is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
// Arguments    : const coder::array<double, 2U> &ibi
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                double AR_order
//                double window
//                double noverlap
//                double nfft
//                double fs
//                struct0_T *output
// Return Type  : void
//
void freqDomainHRV(const coder::array<double, 2U> &ibi, const double VLF[2],
                   const double LF[2], const double HF[2], double AR_order,
                   double window, double noverlap, double nfft, double fs,
                   struct0_T *output)
{
  coder::array<double, 2U> b_z;
  coder::array<double, 1U> F;
  coder::array<double, 1U> b_y;
  coder::array<double, 1U> t;
  coder::array<double, 1U> twt;
  coder::array<double, 1U> x;
  coder::array<double, 1U> y;
  coder::array<double, 1U> z;
  double d_tmp;
  double w;
  double xtmp;
  int i;
  int k;
  int n;
  int nd2;
  int y_tmp;
  if (!isInitialized_freqDomainHRV) {
    freqDomainHRV_initialize();
  }
  // check input
  n = ibi.size(0);
  t.set_size(ibi.size(0));
  for (i = 0; i < n; i++) {
    t[i] = ibi[i];
  }
  // time (s)
  // ibi (s)
  // convert ibi to ms
  // assumes ibi units are seconds
  w = fs / 2.0;
  // prepare y
  n = ibi.size(0);
  y.set_size(ibi.size(0));
  for (i = 0; i < n; i++) {
    y[i] = ibi[i + ibi.size(0)] * 1000.0;
  }
  coder::detrend(y);
  xtmp = coder::blockedSummation(y, y.size(0)) / static_cast<double>(y.size(0));
  n = y.size(0);
  for (i = 0; i < n; i++) {
    y[i] = y[i] - xtmp;
  }
  // Welch FFT
  // calFFT - Calculates the PSD using Welch method.
  //
  // Inputs:
  // Outputs:
  // Prepare y
  d_tmp = 1.0 / fs;
  if (std::isnan(ibi[0]) || std::isnan(d_tmp) ||
      std::isnan(ibi[ibi.size(0) - 1])) {
    output->lomb.f.set_size(1, 1);
    output->lomb.f[0] = rtNaN;
  } else if ((d_tmp == 0.0) ||
             ((ibi[0] < ibi[ibi.size(0) - 1]) && (d_tmp < 0.0)) ||
             ((ibi[ibi.size(0) - 1] < ibi[0]) && (d_tmp > 0.0))) {
    output->lomb.f.set_size(1, 0);
  } else if ((std::isinf(ibi[0]) || std::isinf(ibi[ibi.size(0) - 1])) &&
             (std::isinf(d_tmp) || (ibi[0] == ibi[ibi.size(0) - 1]))) {
    output->lomb.f.set_size(1, 1);
    output->lomb.f[0] = rtNaN;
  } else if (std::isinf(d_tmp)) {
    output->lomb.f.set_size(1, 1);
    output->lomb.f[0] = ibi[0];
  } else if ((std::floor(ibi[0]) == ibi[0]) && (std::floor(d_tmp) == d_tmp)) {
    output->lomb.f.set_size(
        1, static_cast<int>((ibi[ibi.size(0) - 1] - ibi[0]) / d_tmp) + 1);
    n = static_cast<int>((ibi[ibi.size(0) - 1] - ibi[0]) / d_tmp);
    for (i = 0; i <= n; i++) {
      output->lomb.f[i] = ibi[0] + d_tmp * static_cast<double>(i);
    }
  } else {
    coder::eml_float_colon(ibi[0], d_tmp, ibi[ibi.size(0) - 1], output->lomb.f);
  }
  // time values for interp.
  n = ibi.size(0);
  x.set_size(ibi.size(0));
  for (i = 0; i < n; i++) {
    x[i] = ibi[i];
  }
  twt.set_size(output->lomb.f.size(1));
  n = output->lomb.f.size(1);
  for (i = 0; i < n; i++) {
    twt[i] = output->lomb.f[i];
  }
  b_y.set_size(y.size(0));
  n = y.size(0);
  for (i = 0; i < n; i++) {
    b_y[i] = y[i];
  }
  i = ibi.size(0);
  z.set_size(twt.size(0));
  n = twt.size(0);
  for (y_tmp = 0; y_tmp < n; y_tmp++) {
    z[y_tmp] = 0.0;
  }
  if (twt.size(0) != 0) {
    k = 0;
    int exitg1;
    do {
      exitg1 = 0;
      if (k <= i - 1) {
        if (std::isnan(ibi[k])) {
          exitg1 = 1;
        } else {
          k++;
        }
      } else {
        if (ibi[1] < ibi[0]) {
          y_tmp = i >> 1;
          for (nd2 = 0; nd2 < y_tmp; nd2++) {
            xtmp = x[nd2];
            n = (i - nd2) - 1;
            x[nd2] = x[n];
            x[n] = xtmp;
          }
          if (y.size(0) > 1) {
            n = y.size(0) - 1;
            nd2 = y.size(0) >> 1;
            for (k = 0; k < nd2; k++) {
              xtmp = b_y[k];
              y_tmp = n - k;
              b_y[k] = b_y[y_tmp];
              b_y[y_tmp] = xtmp;
            }
          }
        }
        coder::interp1SplineOrPCHIP(b_y, twt, z, x);
        exitg1 = 1;
      }
    } while (exitg1 == 0);
  }
  // cubic spline interpolation
  b_z.set_size(1, z.size(0));
  n = z.size(0);
  for (i = 0; i < n; i++) {
    b_z[i] = z[i];
  }
  xtmp = coder::combineVectorElements(b_z) / static_cast<double>(z.size(0));
  // remove mean
  // Calculate PSD
  b_z.set_size(1, z.size(0));
  n = z.size(0);
  for (i = 0; i < n; i++) {
    b_z[i] = z[i] - xtmp;
  }
  coder::pwelch(b_z, window, noverlap, nfft * 2.0 - 1.0, fs, output->welch.psd,
                F);
  // uses a hamming window
  calcAreas(F, output->welch.psd, VLF, LF, HF, &output->welch.hrv);
  output->welch.f.set_size(F.size(0));
  n = F.size(0);
  for (i = 0; i < n; i++) {
    output->welch.f[i] = F[i];
  }
  // AR
  // calAR - Calculates the PSD using Auto Regression model.
  //
  // Inputs:
  // Outputs:
  // Prepare y
  if (std::isnan(ibi[0]) || std::isnan(d_tmp) ||
      std::isnan(ibi[ibi.size(0) - 1])) {
    output->lomb.f.set_size(1, 1);
    output->lomb.f[0] = rtNaN;
  } else if ((d_tmp == 0.0) ||
             ((ibi[0] < ibi[ibi.size(0) - 1]) && (d_tmp < 0.0)) ||
             ((ibi[ibi.size(0) - 1] < ibi[0]) && (d_tmp > 0.0))) {
    output->lomb.f.set_size(1, 0);
  } else if ((std::isinf(ibi[0]) || std::isinf(ibi[ibi.size(0) - 1])) &&
             (std::isinf(d_tmp) || (ibi[0] == ibi[ibi.size(0) - 1]))) {
    output->lomb.f.set_size(1, 1);
    output->lomb.f[0] = rtNaN;
  } else if (std::isinf(d_tmp)) {
    output->lomb.f.set_size(1, 1);
    output->lomb.f[0] = ibi[0];
  } else if ((std::floor(ibi[0]) == ibi[0]) && (std::floor(d_tmp) == d_tmp)) {
    output->lomb.f.set_size(
        1, static_cast<int>((ibi[ibi.size(0) - 1] - ibi[0]) / d_tmp) + 1);
    n = static_cast<int>((ibi[ibi.size(0) - 1] - ibi[0]) / d_tmp);
    for (i = 0; i <= n; i++) {
      output->lomb.f[i] = ibi[0] + d_tmp * static_cast<double>(i);
    }
  } else {
    coder::eml_float_colon(ibi[0], d_tmp, ibi[ibi.size(0) - 1], output->lomb.f);
  }
  // time values for interp.
  n = ibi.size(0);
  twt.set_size(ibi.size(0));
  for (i = 0; i < n; i++) {
    twt[i] = ibi[i];
  }
  coder::interp1(twt, y, output->lomb.f, output->lomb.psd);
  b_y.set_size(output->lomb.psd.size(1));
  n = output->lomb.psd.size(1);
  for (i = 0; i < n; i++) {
    b_y[i] = output->lomb.psd[i];
  }
  // cubic spline interpolation
  xtmp = coder::blockedSummation(b_y, b_y.size(0)) /
         static_cast<double>(b_y.size(0));
  n = b_y.size(0);
  for (i = 0; i < n; i++) {
    b_y[i] = b_y[i] - xtmp;
  }
  // remove mean
  coder::hamming(static_cast<double>(b_y.size(0)), twt);
  if (b_y.size(0) == twt.size(0)) {
    n = b_y.size(0);
    for (i = 0; i < n; i++) {
      b_y[i] = b_y[i] * twt[i];
    }
  } else {
    b_times(b_y, twt);
  }
  // hamming window
  // Calculate PSD
  // Method 1
  //      [A, variance] = arburg(y,AR_order); %AR using Burg method
  //      [H,F] = freqz(1,A,nfft,fs);
  //      PSD=(abs(H).^2).*(variance/fs); %malik, p.67
  // Method 2
  coder::pburg(b_y, AR_order, nfft * 2.0 - 1.0, fs, twt, F);
  // Method 3
  //       h=spectrum.burg;
  //       hpsd = psd(h, y, 'NFFT', nfft, 'Fs', 2);
  //       F=hpsd.Frequencies;
  //       PSD=hpsd.Data;
  calcAreas(F, twt, VLF, LF, HF, &output->ar.hrv);
  output->ar.f.set_size(F.size(0));
  n = F.size(0);
  for (i = 0; i < n; i++) {
    output->ar.f[i] = F[i];
  }
  output->ar.psd.set_size(twt.size(0));
  n = twt.size(0);
  for (i = 0; i < n; i++) {
    output->ar.psd[i] = twt[i];
  }
  // Lomb
  // calLomb - Calculates the PSD using Lomb-Scargle method.
  //
  // Inputs:
  // Outputs:
  // Calculate PSD
  coder::linspace(w - w / nfft, nfft, output->lomb.f);
  //
  //   [Pn] = lomb(t, y, f)
  //
  //   Uses Lomb's method to compute normalized
  //   periodogram values "Pn" as a function of
  //   supplied vector of frequencies "f" for
  //   input vectors "t" (time) and "y" (observations).
  //   Also returned is probability "Prob" of same
  //   length as Pn (and f) that the null hypothesis
  //   is valid. If f is not supplied it assumes
  //   f =  [1/1024 : 1/1024 : 0.5/min(diff(t))];
  //   x and y must be the same length.
  //  See also:
  //
  //  [1] N.R. Lomb, ``Least-squares frequency analysis of
  //  unequally spaced data,'' Astrophysical and Space Science,
  //  (39) pp. 447--462, 1976.   ... and
  //
  //  [2] J.~D. Scargle, ``Studies in astronomical time series analysis.
  //  II. Statistical aspects of spectral analysis of unevenly spaced data,''
  //  Astrophysical Journal, vol. 263, pp. 835--853, 1982.
  //
  //  [3] T. Thong, "Lomb-Welch Periodogram for Non-uniform Sampling",
  //  Proceedings for the 26th anual international conference of the IEEE EMBS,
  //  Sept 1-5, 2004.
  //  check inputs
  //  subtract mean, compute variance, initialize Pn
  xtmp = coder::blockedSummation(y, y.size(0)) / static_cast<double>(y.size(0));
  z.set_size(y.size(0));
  n = y.size(0);
  for (i = 0; i < n; i++) {
    z[i] = y[i] - xtmp;
  }
  output->lomb.psd.set_size(1, output->lomb.f.size(1));
  n = output->lomb.f.size(1);
  for (i = 0; i < n; i++) {
    output->lomb.psd[i] = 0.0;
  }
  // 	now do main loop for all frequencies
  i = output->lomb.f.size(1);
  for (nd2 = 0; nd2 < i; nd2++) {
    w = 6.2831853071795862 * output->lomb.f[nd2];
    if (w > 0.0) {
      d_tmp = 2.0 * w;
      n = t.size(0);
      twt.set_size(t.size(0));
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        twt[y_tmp] = d_tmp * t[y_tmp];
      }
      x.set_size(twt.size(0));
      n = twt.size(0);
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        x[y_tmp] = twt[y_tmp];
      }
      n = twt.size(0);
      for (k = 0; k < n; k++) {
        x[k] = std::sin(x[k]);
      }
      n = twt.size(0);
      for (k = 0; k < n; k++) {
        twt[k] = std::cos(twt[k]);
      }
      xtmp = rt_atan2d_snf(coder::blockedSummation(x, x.size(0)),
                           coder::blockedSummation(twt, twt.size(0))) /
             2.0 / w;
      n = t.size(0);
      twt.set_size(t.size(0));
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        twt[y_tmp] = w * (t[y_tmp] - xtmp);
      }
      x.set_size(twt.size(0));
      n = twt.size(0);
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        x[y_tmp] = twt[y_tmp];
      }
      n = twt.size(0);
      for (k = 0; k < n; k++) {
        x[k] = std::cos(x[k]);
      }
      F.set_size(twt.size(0));
      n = twt.size(0);
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        F[y_tmp] = twt[y_tmp];
      }
      n = twt.size(0);
      for (k = 0; k < n; k++) {
        F[k] = std::sin(F[k]);
      }
      n = z.size(0);
      if (z.size(0) == x.size(0)) {
        x.set_size(z.size(0));
        for (y_tmp = 0; y_tmp < n; y_tmp++) {
          x[y_tmp] = z[y_tmp] * x[y_tmp];
        }
      } else {
        times(x, z);
      }
      d_tmp = coder::blockedSummation(x, x.size(0));
      n = z.size(0);
      if (z.size(0) == F.size(0)) {
        F.set_size(z.size(0));
        for (y_tmp = 0; y_tmp < n; y_tmp++) {
          F[y_tmp] = z[y_tmp] * F[y_tmp];
        }
      } else {
        times(F, z);
      }
      w = coder::blockedSummation(F, F.size(0));
      x.set_size(twt.size(0));
      n = twt.size(0);
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        x[y_tmp] = twt[y_tmp];
      }
      n = twt.size(0);
      for (k = 0; k < n; k++) {
        x[k] = std::cos(x[k]);
      }
      n = twt.size(0);
      for (k = 0; k < n; k++) {
        twt[k] = std::sin(twt[k]);
      }
      n = x.size(0);
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        xtmp = x[y_tmp];
        x[y_tmp] = xtmp * xtmp;
      }
      F.set_size(twt.size(0));
      n = twt.size(0);
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        xtmp = twt[y_tmp];
        F[y_tmp] = xtmp * xtmp;
      }
      output->lomb.psd[nd2] =
          d_tmp * d_tmp / coder::blockedSummation(x, x.size(0)) +
          w * w / coder::blockedSummation(F, F.size(0));
    } else {
      n = z.size(0);
      if (z.size(0) == t.size(0)) {
        x.set_size(z.size(0));
        for (y_tmp = 0; y_tmp < n; y_tmp++) {
          x[y_tmp] = z[y_tmp] * t[y_tmp];
        }
      } else {
        times(x, z, t);
      }
      d_tmp = coder::blockedSummation(x, x.size(0));
      n = t.size(0);
      x.set_size(t.size(0));
      for (y_tmp = 0; y_tmp < n; y_tmp++) {
        xtmp = t[y_tmp];
        x[y_tmp] = xtmp * xtmp;
      }
      output->lomb.psd[nd2] =
          d_tmp * d_tmp / coder::blockedSummation(x, x.size(0));
    }
  }
  //  return denormalized spectrum (see T. Thong)
  output->lomb.psd.set_size(1, output->lomb.psd.size(1));
  n = output->lomb.psd.size(1) - 1;
  for (i = 0; i <= n; i++) {
    output->lomb.psd[i] = output->lomb.psd[i] / static_cast<double>(y.size(0));
  }
  // calc lomb psd
  b_z.set_size(1, output->lomb.psd.size(1));
  n = output->lomb.psd.size(0) * output->lomb.psd.size(1) - 1;
  for (i = 0; i <= n; i++) {
    b_z[i] = output->lomb.psd[i];
  }
  calcAreas(output->lomb.f, b_z, VLF, LF, HF, &output->lomb.hrv);
}

//
// File trailer for freqDomainHRV.cpp
//
// [EOF]
//
